Illusions in the Ring Model of visual orientation selectivity 
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Abstract 

The Ring Model of orientation tuning is a dynamical model of a hypercolumn of visual area VI in 

the human neocortex that has been designed to account for the experimentally observed orientation 

tuning curves by local, i.e., cortico-cortical computations. The tuning curves are stationary, i.e. time 

independent, solutions of this dynamical model. One important assumption underlying the Ring Model is 

that the LGN input to VI is weakly tuned to the retinal orientation and that it is the local computations 

in VI that sharpen this tuning. Because the equations that describe the Ring Model have built-in 

l/~) equivariance properties in the synaptic weight distribution with respect to a particular group acting on 

the retinal orientation of the stimulus, the model in effect encodes an infinite number of tuning curves 

that are arbitrarily translated with respect to each other. By using the Orbit Space Reduction technique 

Xy^ we rewrite the model equations in canonical form as functions of polynomials that are invariant with 

^\ respect to the action of this group. This allows us to combine equivariant bifurcation theory with an 

efficient numerical continuation method in order to compute the tuning curves predicted by the Ring 

Model. Surprisingly some of these tuning curves are not tuned to the stimulus. We interpret them as 

neural illusions and show numerically how they can be induced by simple dynamical stimuli. These neural 

illusions are important biological predictions of the model. If they could be observed experimentally this 

would be a strong point in favor ot the Ring Model. We also show how our theoretical analysis allows 

to very simply specify the ranges of the model parameters by comparing the model predictions with 

£> published experimental observations. 

m 

On 

Tt Author Summary 

r-»«. The visual cortex of humans features a columnar organization. Each column contains cells whose receptive 

fields monitor almost identical retinal positions and orientations. Such nearby orientation columns encode 
different orientations. A set of orientation columns encoding all possible orientations is a hypercolumn of 
orientation. The visual input to a hypercolumn is from the retina through the thalamus. The hypercolumn 
^- is a relay that sharpens the output of the thalamus to make the local retinal visual orientation more 

" conspicuous" . This happens because the neurons of the hypercolumn inhibit and excite each other 
according to simple rules. These rules are used to build a mathematical model of a hypercolumn whose 
predictions can be compared with experimental evidence. We study the properties of such a model whose 
rich symmetries can be traced to the fact that two orientations differing by the value it are the same 
visually. Important consequences of our analysis are the clarification of the model parameters role in 
shaping the "perception" of the hypercolumn and the prediction of a neuronal illusion, the fact that a 
hypercolumn can be in the same state as if the retinal stimulus were at an orientation rotated by 90 
degrees with respect to the actual one. 

Introduction 

Since the discovery by Hubel and Wiesel [1] of the selective response of a single neuron to some orien- 
tations, a long standing debate has been the degree of cortical computation involved in this selectivity 
compared to the feedforward selectivity implied by the LGN projections. Cortical models [2, 3, 4] have 
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Figure 1. A simplified view of the visual path from the retina through the LGN to cortical area VI. 
The receptive field of the LGN cells feeding the hypercolumn of orientation of VI contains a grating of 
orientation xq. This grating excites mostly the LGN cells that share this receptive field and are aligned 
in the direction xo; the tuning is broad, see the curve in the middle part of the figure. These LGN cells 
project onto the network of cells in the hypercolumn of orientation of VI whose interactions, 
represented by the Ring Model, result in a sharpening of the tuning around the grating direction, see 
the curve in the righthand side of the figure. 



been used to show how this selectivity can be produced in a cortex with center-surround interactions in 
the orientation domain. 

The Ring Model of orientation tuning was introduced by Hansel and Sompolinski [4] and studied by 
several other scientists [5, 6, 7, 8, 9], after the seminal work of Ben-Yishai and colleagues [3], as a model of 
a hypercolumn in primary visual cortex. This rate model is a simplification of complex spiking networks 
[2, 10] designed to make it easier to understand the role of different mesoscopic parameters. 

It assumes that the local orientation x in the receptor fields of the neurons in the column is encoded in 
their activity, or firing rate, noted A{x, i). The interaction between the neurons is modeled by a function 
J of the orientation that represents how the activities corresponding to two different orientations reinforce 
or inhibit each other. This function is called the connectivity function of the model. With this in mind, 
the dynamics of the firing rate can be represented by the following integro-differential equation in the 
line of the model of Wilson-Cowan [11]: 



rA(x,t) = -A(x,t) + S 
A(x,0) = A (x) 



tt/2 

A J J(x-y)A(y,t))dy/n + eI(x)- 

V -tt/2 



t > 



(1) 



t defines the intrinsic dynamics of the population and / is the input from the lateral geniculate nucleus 
(LGN) to the hypercolumn whose contrast is defined by the parameter e, see figure 1. S is the sigmoid 



S{x) 
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1 



which takes values between and 1, A is a parameter that determines the nonlinear gain of the sigmoid, 
is a threshold that controls for which value the sigmoid takes the value 1/2, the function / represents 
the input from the LGN. 

Regarding the connectivity function, in all cases it is an even 7r-periodic function, positive for ori- 
entation values close to (corresponding to an excitation) and negative for orientation values larger in 
magnitude (corresponding to an inhibition). The rich symmetries of J play a prominent role in our 
upcoming analysis of the Ring Model. The reason for this is that, when the contrast e is equal to 0, 



equation (1), is equivariant, i.e. it has some nice properties with respect to the action of a certain group 

that we proceed to describe. 

Let us consider the group, noted T, of translations of ~R/2irZ. An element T 7 , 7 € M/27rZ of this 

group acts on the orientation x by T 7 • x = 7/2 + x and on the activity function A(x, t) by T 7 • -A(x, t) = 

A(7/2+£, t). Similarly, we consider the reflection, noted R, such that R-x = —x and therefore R-A(x, t) = 

A{— x, t). We note G the group generated by T and R. We can abstractly rewrite equation (1) as T(x) = 0, 

[ / *7 2 " 

where T(x) = rA(x,t) + A(x,t) - S A( J J(x - y)A(y,t))dy/-w + el{x) - 6) 

V -7r/2 ' _ 

properties of the function J it is easy to verify that T satisfies the condition T{g ■ x) — g ■ J-(x) for all 
elements g in the group G when e = 0. Since moreover we clearly have 



Using the symmetry 



T^T J2 = T J1+r , RT^ = T_ 7 i? for all 71, 72, 7 S E./2irZ 
T = Id R 2 = Id 

the group G is in effect isomorphic to 0(2), the group of two-dimensional othogonal transformations [12, 
Chapter 1]. 

Several variants of this model have been studied in the literature, e.g., in [8, 9] J is a difference of 
Gaussians while in [•>] the authors start with a network of excitatory/inhibitory spiking neurons and 
derive a meanfield approximation of this network yielding an interaction function J described by the 
following equation: 

J(x) = J + Ji cos(2x) (2) 

The input I from the LGN has a similar shape 

I(x) = 1-13 + 13 cos(2(a; - x )). (3) 

As mentioned above, and as shown in figure 1, it is weakly tuned, i.e., maximal, at x = xq and it is 
the network, modeled by (i), that sharpens this tuning. These authors vary the anisotropy parameter j3 
between and 1. 

The sigmoid S is often chosen to be a Heaviside function, or, as in [ ], a piecewise linear approximation 
of the sigmoid. In [6, 7, 8] it is a true sigmoidal function. 

The parameter J\ is positive, an important property of the network that is necessary to produce the 
tuning curves.. Jo is most of the time negative [3, 7, 8, 9] but can be positive as well [9]. Note that the 
JiS, i = 0,l are the first Fourier coefficients of the function J. Jo is its mean value and can be positive 
even if the surround is inhibitory. For example, in [ ], we find Jo = —7.3, J\ = 11, j3 = 0.1, 9 = 
which are the values in [3] except for 9 = 1 and e = 0.164. The slope, or nonlinear gain, A is assumed 
to be equal to 1. Using the previous rescaling, it becomes J = —1, J\ = 1.5, A = 7.3/S"(0) = 29.2 and 
9 — > 9/7.3 « 0.1 in the case of [3] and 9 = in the case of [ ]. 

The model (1) is called an activity model in the terminology of [()]. For technical reasons we turn 
it into a voltage model as follows. We first rewrite equation (1) in a more compact and convenient, 
functional, form : 

tA = -A + S [A( J -A + eI-9)}. 

J is now thought of as a linear (integral) operator acting on the function A as the periodic convolution 

tt/2 

J • A(x,t) = J J(x — y)A(y,t))dy/ir, see, e.g., [13]. We then perform the change of variable V = 

-tt/2 

J ■ A + el — 9. Assuming that the input current is not a function of time, this leads to the following 
equation 

tV = ~V + J-S(\V)+eI-9 (4) 

Note that this equation, as (1), is G-equivariant. 



The stationary solutions (some of them called tuning curves, see the paragraph Discussion) of (1) 
(respectively of (4)) satisfy ^4 = (respectively V = 0). Characterizing and computing them for different 
values of the parameters is the first step toward understanding the dynamics of the solutions to these 
equations. Indeed, it is known that this type of equations is such that there are only heteroclinic (linking 
two stationary solutions, or equilibria) or unbounded orbits. Since we showed in [ "] that all trajectories 
were bounded, this implies that they are made of heteroclinic orbits. This motivates further the study 
of the stationary solutions of (4) . One of our goals is to show how the stationary solutions are organized 
and to give indications about the dynamics in a given range of parameters, corresponding to biologically 
plausible values. This is relevant because many large scale models of VI including many hypcrcolumns 
represent them with the Ring Model. Therefore a good understanding of one hypercolumn paves the way 
to an understanding of a population thereof. We show that, depending on the nonlinear gain A, there 
may exist many stationary solutions, which are all acceptable responses of the network to a given input 
from the LGN (at least for the model at hand). Thus this local orientation tuning device may behave 
less trivially than what it was initially designed for. In effect, the existence of these stationary solutions, 
sometimes called persistent cortical states, can make the local dynamics quite intricate when A is large 
enough to support the existence of these extra solutions. 

It is worth noticing that the stationary solutions of (1) and (4) are in one to one correspondence. As a 
consequence we will work on (4) because it is mathematically more convenient. 

We will follow a method similar to the one developed in [ ; I] to compute the stationary states of (4) . 
The method has been modified to take into account the symmetries of the Ring Model. The general 
idea is that the LGN input is weak and only modulates the network activity. Hence the cortical network 
(represented by the Ring Model) encodes the possible tuning curves within its connectivity function and 
when presented with a weak external input, produces small deviations of 'its' tuning curves. Our goal 
is to compute these tuning curves. However, because of the symmetries of the connectivity function J, 
the model in effect encodes an infinite number of tuning curves, and this is an endless cause of numerical 
problems. Indeed we pointed out above that if the input current was null, equation (4) (respectively 
(I)) was G-equi variant. This implies that if V(x) is a stationary solution of (4), so are V(x + j/2), 
7 £ ~R/2irZ and V(—x). We show that by performing an appropriate change of variables, we can get rid 
of this redundancy and recover numerical accuracy. 

Materials and Methods 

Turning the problem into a finite dimensional one 

Problem (1) (respectively (4)) is infinite dimensional since the solutions live in some (unspecified, but 
a priori infinite dimensional) functional space. By truncating the Fourier series of the even 7r-periodic 
connectivity function J we reduce the problem to a finite number of dimensions. We write 

JV 

J(x) — J + y_] J P cos(2px) 
P =i 

where N is the number of Fourier modes that are sufficient to represent J. We will show how the choice of 
N affects the biological functional properties of the Ring Model. Notice that by varying N, we generate 
a family of models that contains all previously published ones. 

Remark 1 For convenience, we shall write cosfc for the junction x — > cos(fca;). The same holds for sin^. 

It was shown in [13] that this form of the connectivity function implies that the solutions to (4) can be 
written V(x,t) = V^(x,t) + V ± (x,t), where V^ is a linear combination of the functions C0S2 P and ski.2p, 
p = 0, ••• ,N and the function V 1 - tends to exponentially fast when t-> oo. This implies that the 
stationary solutions satisfy V 1 - = 0. 



Remark 2 The spectrum of the integral operator J associated with the eponymic connectivity function is 
readily seen to be equal to S(J) = < Jq, -# f 

I ■* J l<p<N 

We can, up to a rescaling of A in (1), assume that Jo takes the values ±1: 

J = £ e {-1,1}. 

Similarly we define £*, = ±1 by 

Jk = e k \J k \,e k e{-l,l} k = l,---,N 
With all this in hands, the tuning curves satisfy the equation: 



tt/2 



V(x) 



-7T/2 



/V 



Jo + ^J P cos 2p (x - y) 



S[\V(y,t))]dy/ir + eI(x) 



(5) 



It follows that any stationary solution to (4) can be written 



N 



V(x) = v a + ^2 \J\ J p\ v v^ cos 2p (a;) + v { p 2) sin 2p (a;) 
P =i 



,(1),(2) 



where vq,v p , p — 1,- ■ ■ ,N are 2N + 1 reals. Solving (5) is therefore equivalent to finding these reals. 
In the case of a general solution, V"(x,t) is given by the same formula where the coefficients are now 
real functions of time. Under the assumption that V^ is neglected, it is easy to obtain the system of 



ordinary differential equations that are satisfied by the functions Vq,v p 



(i),(2) 



Using the complex values 



Zk 



dcf (i) 



+ ivl z >, fc = l, 



vo + vo 



, N we obtain the following equations 



f / S 



N 



A«o + A E ^PR (z p e~ 2 P*y) 
P =i 



dy — 9 + eIq 



dcf 



= B (v ,{zp}) -0 + sIq 

N 



Zk + Zk = £k- 



'\Jk\ 



IS 



Awo + A £ V\Jp\n {z p t 



-2piy\ 



p=l 



(6) 



a 2kiy 



dy + elk 



dcf 



= -B fc (^o, {> P }) + elk k = 1, • 



,iV 



dcf 



JV 



where /(«) = 7 + E Wl J fcl e and I el, I* e 



fc=i 



The coefficients defining the tuning curves satisfy the following equations: 



v Q = B (v , {z p }) -6 + el 

Zk = B k (v ,{z p }) +el k k = 1, 



,7V 



The N + 1 dimensional vector (v , Z\,- ■ ■ , Zjsr) is a representation of V" . The group G also acts on this 
representation as follows 



T 7 • (v ,zi,z 2 ,--- ,z N ) 
R- (v ,z 1 ,z 2 ,- ■■ ,z N ) 



(v ,e 2 ^ Zl ,e^z 2l --- ,e 2m ~<z N ) 7 G 

(v ,Zl,Z2," ■ iZn) 



As shown in the introduction, the use of the group 0(2) is motivated by the fact that when e = 0, then 
(6) are 0(2)-equivariant 1 . Since if V? is a stationary solution of (6) for e = 0, so is g ■ V* , V<? G 0(2) 



x If we write (6) when e = as 4vll = F(VH), this means that Vg e 0(2),W', F(g ■ V"H) = g-F(Vll) 



N 

there is an infinity of tuning curves that are encoded by the network. However, when £ Y\ Ik ^ 0, all the 

fe=i 
symmetries are broken, (6) are not 0(2)-equivariant anymore and the number of tuning curves becomes 
finite. 

In the next two sections we study the cases N = 1 and N = 2. The second case shows that adding 
more modes does not change the main results of the analysis. 



Keeping only one mode in J 

We consider the following connectivity function 

J = £ + J\ cos 2 , J\ > 

From our previous analysis of the symmetries of the Ring Model, we know that the equations (6) are 
redundant when e = 0. In order to eliminate this redundancy they should be rewritten using their 
equivariant structure with respect to the action of the group 0(2). In this case, this turns out to be 
equivalent to writing an equation for vq and the magnitude p of zi. It is convenient to write z\ = 
Vi + iv2 = pe 2l,p , which yields to the following equations, assuming xq = 0: 



v = -vq + e B {v , p) -6 + e(l - j3) 
P = -p + JiB 1 (vo,p) + ^=cos 2 {(p) 
2pip = -sin 2 (<^)^§= 



(7) 



for the dynamics, and 



+ e(l-/3) 
el 
V7I 



eP cos 2 (^) 



(8) 



V = £ Q B (v ,p)- 

P = VJiB 1 (v ,p) 
= OnaMft 

for the tuning curves. The functions B and B\ are given by (using an integration by parts to factor out 
pin Bi{v ,p)) 

B (vo,p) = £ ] S(X(v Q + ^hp cos 2 x))f 



Bi(v ,p) 



dx 



rj\\p J S'(\(v + \f7\p C0S2X)) sin 2 a; ^ 



Equations (7) do not produce the same dynamics as (6) because the change from Cartesian to polar 
cordinatcs is not a diffeomorphism. Nevertheless equations (8) are most useful for computing the tuning 
curves. 



Case of 2 modes, N = 2 

We now write: 



J = Eq + J i COS 2 + J 2 COS4 



In order to agree with known experimental facts, the tuning curves should be mainly unimodal. Compared 
to the previous case, the fact that the second mode is nonzero could induce an "interaction" between the 
two modes leading to multimodal tuning curves. 

Following the analysis of section we have to solve five coupled equations which are redundant because 
of the action of the group G which in this case reads 



T 7 - (V , Z!,Z 2 ) = 
R- (uo,Zl,Z2) = 



K,e 2 ^z 1 ,e 4 ^z 2 ) 
(v ,zi,z 2 ) 



In order to eliminate the redundancy arising from this symmetry we used polar coordinates as in the case 
N = 1. It is tempting to do the same with the two complex variables Z\ and z 2 but it turns out to be a 
dead end. 

The main reason is numerical: we compute (see Supporting information, second paragraph) the solu- 
tions of the nonlinear equations (6) using numerical continuation. This scheme works well if the Jacobian 
of the nonlinear equation has at worst a one-dimensional kernel at some isolated points. If we write 
z i = Pie 2lVl , z 2 = /3 2 e 4lV2 , then, using the invariance by T 7 , the equations for ipi = 0, i = 1, 2 simplify to 
= Fi(i>o, pxipii^px — ^2) which are functions of the phase difference y>i—ip2- It turns out that the other 
equations (for va,p\,p2) also involve only tpi — if2- We were unable to find a simple relation between 
.Fi and F%. As a consequnce we end up with 5 equations in the 4 unknows Vq, pi, p2,fi — f2'- this is 
unappropriate for numerical continuation and we need to find a way to obtain 4 equations in 4 unknowns. 

To reach this goal we turn to a general technique, the Orbit Space Reduction [14], which provides the 
right change of coordinates through the use of what is called a Hilbert Basis 2 . A fundamental property 
is that any smooth equivariant function can be expressed using the elements of a Hilbert Basis and their 
gradients. 

A Hilbert basis associated to the action of the group 0(2) ([14, page 205]) is given by the 0(2)-invariant 
polynomials: 

7ri = Z1Z1, ir 2 = z 2 z 2 , 7T 3 = 3? (zlz 2 ) 

which must satisfy the constraints 

7T1 > 0, 7T 2 > 0, 7T 2 < 7T 2 7r 2 (9) 

Our analysis is now focused on the so-called Orbit Space, i.e. the subset of K 4 of the four-tuples (vq,7t), 
7? = (tti, 7T2, 7r 3 ), that satisfy the previous inequalities. 

As the fonction Bq(vq, z\, z 2 ) is (9(2)-invariant (i.e. Bo(g ■ (vq, z\, z 2 )) = Bq(vq, z\, z 2 ) for all g in G), 
it is a function, noted Bq(vo, n), of the variables (vq, ff), i.e. Bq(vq, z\, z 2 ) = Bq(vq, n). Furthermore since 
the pair (Bi,B 2 ) is 0(2)-equi variant (i.e. g ■ (B x , B 2 ){v , z%, z 2 ) = {B x {g ■ (v , Z\, z 2 )), B 2 (g ■ (v ,z 1 ,z 2 ))) 
for all g in G), it can be written: 

-Z\ +B 1 (v ,z 1 ,z 2 ) =a(v ,Tr)z 1 + b(v , Tr)z 1 z 2 , . 

-z 2 + B 2 (v ,z 1 ,z 2 ) = c(v ,tt)z 2 + d(v ,Tf)zf 

where a(vo, n), b(vo, ff), c(vq, tt), oI(vq, if) are 0(2)-invariant functions. Notice that this implies that Bi(vq, 0, z 2 ) 
0. 

Using the definition of the polynomials 7Tj, it is possible to rewrite (6) only in terms of (vo,tt) (as we 
did in the previous section with the polar coordinates): 

"0 = -Vq + B Q (vq,Tt) -9 

Tr 1 = 2a(v ,Tr)Tr 1 +2b(v ,ir)ir 3 ^ 

tt 2 = 2c(w , 7r)vr 2 + 2d(v , 7?)7r 3 

7r 3 = [2a(v , 7?) + c{v , 7?)] 7T 3 + 2c(v Q , 7r)7ri7r 2 + d(v , 7t)7T 2 

These equations are solved to find the tuning curves. 

2 The ring 'R.q of 0(2)-invariant polynomials is finitely generated as an M-algebra, this goes back to Hilbert. A family 
7Ti , ■ ■ ■ , 7r s of generators of TIq is called a Hilbert Basis. 



Results 

We used the equations of the previous paragraph to find the tuning curves in the cases N = 1 and N = 2. 
Some of these tuning curves corresponded to a neuronal illusion. We then showed that adding more 
modes to the connectivity function did not change the results. Finally we designed two different types of 
external stimuli for bringing the network to the illusory states. 

Finding the tuning curves, case N = I 

The Ring Model is based on one main ingredient: at null contrast and for small values of the nonlinear 
gain A, there is a unique stationary solution, which is not tuned! Indeed, this stationary solution has to 
satisfy p? — otherwise it would not be unique because of the group 0{2) equivariance. Thus, in order 
to produce tuning curves (that are tuned by definition), we need a solution to (8) satisfying pf ^ 0. This 
means that we must investigate for which values of A, if any, the p solution of (8) bifurcates. For no 
external input they read: 

pf = y/^Biivo,^) (12) 



A bifurcated solution arises when 



Ji f <-.//, / « rr \ \ . • 2 * Ji 



J 1 d p B 1 (v [) ,p)\ p=0 ^\— I S'(\(v o + 0- V JiCos 2 ))sin 2 = A — S'(Xv ). 

It follows that the equations for the existence of a tuned stationary solution, satisfying p* ^ are: 

4 = e S(\v£)-" 



1 



ASW)4 (13) 



Using the relation S' = S(l — S) it is straightforward to show that A and J\ must satisfy the condition 
AJi > 8. In (see Supporting information, third section), we find other equivalent conditions that are 
used to produce the graphs shown in figure 2. These graphs show that the threshold 9 and the ratio 
excitation/inhibition are constrained in order to produce the tuning curves. When these conditions are 
satisfied we obtain a continuum of tuning curves parametrized by the phase angle </?, noted TC V , which 
are given by 

TC v {x) = S \(d + \y/Jip f cos 2 (2: - <p) - 0) 

Note that these tuning curves are dynamically stable because they are produced by a pitchfork bifurca- 
tion, as can be seen by examing equation (12). The bifurcated branch of interest is the one corresponding 
to p f > 0. 

The next step is to investigate what happens when we switch on the LGN drive, i.e. when e ^ 0. First, 
when the anisotropy j3 of the LGN input is not zero, the equations (6) are not 0(2)-equivariant anymore. 
This is a symmetry breaking and, as mentioned above, there are a finite number of tuning curves. Two 
important questions are 1) how many of the (continuum of) tuning curves remain solutions and 2) what 
is their stability? For very small s ^ 0, switching on the LGN can be viewed as a perturbation of the 
nonlinear equations when e = 0, as a consequence, we expect an opening of the pitchfork as we explained 
in [13]. This is confirmed by figure 3. 

We know from our previous analysis that these solutions satisfy: 
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Figure 2. Boundaries in the (J\, 9) plane for the existence of the tuning curves for e = — 1 (Green, it 
is very close to J\ = 108 + 1), and £q = 1 (Red). The domain of existence lies above these boundaries. 



v{ ) =e B (v f 0: pf)+s(l-fJ)-8 
2cpf = kn, keZ 



Considering the two cases k even and k odd we obtain: 



J - 



e B (v£,pf) + e(l-P)-8 + 



9: 



f 



VJiB^vIpD 

kn, fcgZ 



ef3 



or 



^ = £ -B (^,P^) + £(l-/3) 
J lBl (v f Q ,pf)-^ 



Pi- 

<P f o 



In 
2 



(2k- 



i)f, fcez 



Because p — > B (v ,p) is even and p —> Bi(v Q ,p) is odd, necessarily p{ = —pi- We solve these equations 
for (dq, p e ) as functions of A by using a continuation algorithm described in [13], the results are shown 
in figure 3. 

From this bifurcation diagram, we observe that there are three stationary solutions for A > 10, one 
unstable corresponding to a small value of p (thus it is untuned and, by definition, does not represent a 
tuning curve) and two others which arc tuning curves. One, noted TCq, is peaked at x = and the other, 
noted TC^i2i is peaked at x = tt/2. Note that the values A > 10 are in agreement with the previously 
derived necessary condition AJi > 8. The two tuning curves are shown in figure 4. The interesting fact 
to notice is that the tuning curve TC^/2 is a somewhat bizarre stable state of the Ring Model. We may 
want to call it a neuronal illusion, the 90 degrees illusion, since it corresponds to the fact that, even if the 
thalamic input is peaked at the zero degree orientation, the Ring Model may be (and stay) in a stable 
state corresponding to a tuning curve peaking at 90 degrees! In other words, even if the thalamic input 
"says" degrees, the hypercolumn of orientation "says" 90 degrees. 



Finding the tuning curves, case N = 2 

The previous results may seem to depend very much on the type of simple connectivity function that we 
have assumed. In fact this is not so. By adding one more mode to this function we show that they are 
generic if the resulting function preserves the structure of the local excitation and the lateral inhibition. 
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Figure 3. Plot of (vq, p e , p ) shown in red, green, blue, respectively, as functions of A for 
e = 0.01, 9 = 0, J\ = 1.5, j3 = 0.1. Notice the turning points labelled with black dots. 




Figure 4. The tuning curves TCq and TC^/2 for A = 15 and 
e = 0.01, x Q = 0, 9 = 0, ,h = 1.5, /3 = 0.1. 
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The tuning curves are now solution of the nonlinear equations: 

v = B (v ,tt) -9 

= 2a(v , 7r)7ri + 26(w , 7r)7r 3 

= 2c(v , T?)ir 2 + 2d(v , 7?)7r 3 

= [2a(vo, 7?) + c(v , 7?)] 7r 3 + 2c(v , 7r)7ri7r 2 + d(v , tx)t:\ 

We hrst look at the case where the sigmoid function S is zero at the origin, i.e. Sq(x) — S(x) — h. 



(14) 
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Figure 5. Left: Connectivity function used in the example described in the text. Right: Plot of the 
solutions obtained with Sq instead of S as functions of the nonlinear gain A. Each solution is made of a 
4-tuple (red, green, blue, violet). The parameter values are J = — 1, J\ = 9, J 2 = 6.66, 9 = 0, e = 0. The 
'wavy' branches are numerical artefacts of the approximation of the centered sigmoid Sq by polynomials. 
Note that, as expected because we are on the Orbit Space, the values of 7Ti and tt 2 are positive, see text. 

We see on the graph of solutions on the Orbit Space shown in figure 5 that there are two bifurcations 3 
from the trivial solution (vq,tt) = 0, at the points, noted Pi and P 2 in this figure, corresponding to the 
values Ai < X 2 of the nonlinear gain. Considering again figure 5 it motivates the following remarks: 

1. The first bifurcated branch from Pi reaches high values well before P 2 . 

2. Our Orbit Space reduction procedure allows us to compute numerically such secondary bifurcation 
points as P3 which might produce linearly stable tuning curves. These are undesirable from a 
biological viewpoint because they produce stable multimodal tuning curves. 

The linear stability analysis shows that the branch bifurcating from P\ is stable and corresponds to 
a continuum of tuning curves parametrized by the phase qngle if and given by: 

V^ TC 1 (x) = S 



\(vq+ y tt[ J 1 cos 2 (ie + ip) - 6 j 



The unstable tuning curves bifurcating from P 2 (before P3) are given by: 



A f Vq + yir 2 \J 2 \ cos 4 (x + ip) 



y V TC 2 Jx) = S 



This shows that in order to have unimodal tuning curves (responses of the hypercolumn represented by 
the Ring Model), it is necessary that Ai < X 2 or equivalently J 2 < Ji because the nonlinear gains which 



3 These are not regular bifurcations because we are working on the Orbit Space. 
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produce the pitchforks are «Sg(0)Ai = j-, i = 1,2. Moreover, since tt{ quickly reaches high values, TC\ (0) 
is close to 1: the response does not depend upon the contrast e of the LGN. This implies that the working 
range of the nonlinear gain X is close to the value X\. 

It is now possible to understand the diagram of solutions shown in figure 6 obtained with the regular 
sigmoid 5 as a deformation of the diagram shown in figure 5. As in the case N = 1, the bifurcated 
branches will persist if the coefficients Ji, i = 1,2, and the threshold 9, satisfy the constraints shown in 
figure 2. Indeed, if we reproduce the analysis in the paragraph Materials and Methods, the existence of 
a pitchfork for the Zi, i = 1,2 coordinate is given by 

f v f Q = e S(Xv f )-6 
\ 1 = XS'(Xv f )% 

We again notice that the first branch bifurcating from Pj (in green in figure 6) is quickly reaching 
high values and that the tuning curve is now asymmetric (this is much easier to see in the middle part 
of figure 7) . This is because the ~ki , 773 components (in blue and magenta in figure 6) are not zero as in 
the case N — 1. The tuning curve corresponding to the first bifurcated branch is given 



TC(x) = S 



( 



+ y n{ Ji cos 2 (x) + y^i\J 2 \ cos 4 (x + (p 2 - <pi) 



(15) 



where z\ — jTt\e 



,2iip x 



Z-2 



ir 2 e 



iiip 2 



and ^1,^/71^0084(932 — fi) — i"3- Remember that there is an infinity 



of tuning curves that are obtained from the one given by equation (15) by applying an element of the 
group G. 




Figure 6. Plot of the solutions with the regular sigmoid S as a function of the nonlinear gain A. Case 
J = -1, h = 9, J 2 = 6.66, 6 = e = 0. 



We have plotted in figure 7 examples of the tuning curves for three values of the nonlinear gain A that 
are slightly larger than the values corresponding to the three bifurcation points Pi, P 2 and P3 in figure 
6. These tuning curves are obtained by reading from figure 6 the 4-tuple (vo,tt). This yields, through 
the relation iri^/ir^ cos4(ip 2 — ipi) — ir^, the value of ip 2 — ipi that is needed in equation (15). Notice 
the unstable multimodal tuning curves that appear once the stable tuning curve has saturated (middle 
plot in figure 7). This is an indication that the nonlinear gain should not be too high otherwise most 
responses of the network will be saturated. 
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Figure 7. Tuning curves at slopes A = 0.26 (left), A = 0.355 (middle), and A = 0.58 (right) when the 
input is equal to 0, see text. On the left and in the middle, stable tuning curves are shown in 
continuous line, unstable ones in dotted lines. Stability is not shown in the plot on the right, except for 
the null solution. The other parameters are the same as in figure 6. 

If we switch on the LGN, the contrast e is nonzero. The external current is given by I(x) = Iq + 
JiV ' J\ COS2 +Ii\l\Jz\ COS4. If l\Ii ^ 0, the symmetries of the equation (11) are broken and we expect a 
finite number of solutions. More precisely, the same argument as for N = 1 shows that we can interpret 
(6) as a perturbation of (11) when e is small, leading to an opening of the pitchforks. Hence, when the 
nonlinear gain A is close to that of P\ we have z 2 ~ for e small. Since the equations 

-i>o + B Q (v (h zi,z 2 ) - + sIq = 
-Z\ + Bi(v Q ,zi,z 2 ) +eh = 

are the same as in the case N = 1 when z 2 = 0, they do not change much when z 2 ~ and the 90 degrees 
illusion found in the previous section remains: there are two tuning curves, one peaking as the external 
input I and one translated by 90 degrees. This analysis is confirmed by the results of the numerical 
computations shown in figure 8 where we show the solutions of (6) for N = 2, Iq = 1 — /3, \[J\I\ = 
P, VW\h = 0.1/3. 




Figure 8. The tuning curves TC and TC^/2 for A = 0.26, 

e = 0.01, x = 0, = 0, Ji = 9, J 2 = 6.66, $ = 0.05 and I o = l-0, y/Jlh = P, 
plotted in red. Notice the unstable weakly tuned tuning curve shown in black. 



0.1/3. I is 
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Arbitrary number of modes 

We can perform the same computations using more modes, this will only bring in more tuning curves. 
Because these tuning curves only appear once the stable unimodal tuning curve has saturated these high 
values for the nonlinear gain A are biologically irrelevant or nonplausible. Notice also that the neuronal 
illusions found in the case N = 1, J\ > are still present for N > 1, J\ > 0, as shown for example in 
figures 7 and 8. Indeed, as seen in the previous section, they only depend upon the fact that the network 
features a pitchfork bifurcation at the point noted P\ in figures 5 and 6 and this is always the case for 
any value of the number N of modes if the coefficients Ji satisfy the mild constraints we have described 
previously and we summarize in the next section. 

Dynamical 90 degrees illusions 

In the last paragraphs, we found two cortical representations of the same external stimulus. An immediate 
question is how can we bring a hypercolumn of orientation into each of these two states? Can we drive 
the cortical state to the illusion using only the stimulus I? We answer this question positively in the next 
two paragraphs. 

Rotating the stimulus back and forth 

As the illusory tuning curve TC^/2 is very close to the cortical state corresponding to the response of the 
network to a stimulus peaked at 7r/2, we first present a stimulus peaked at to put the system in the 
TCo state corresponding to no illusion. We then slowly change the position of the peak of the external 
stimulus / and bring it to the value n/2. The network follows the stimulus and its response is peaked 
at 7r/2. We then suddenly change the stimulus to a stimulus peaked at and since the responses of 
the network to stimulus oriented at or ir/2 are very close, the cortical state will remain in the state 
peaked at 7r/2, the one it is in just before the sudden change of the stimulus. This is reminiscent of the 
after-effect illusion and can be confirmed by numerical simulation. 

The resulting effect is shown in figure 9 when the time variation of the peak xq of the stimulus in 
equation (3) is given by 

fmin( T ^,l) iff €[0,2000] 
x (t)={ f if t € [2000, 2e4] 

if t > 2e4 

Using a mixture of the two stimuli 

This second dynamical 90 degrees illusion is very close in principle to the one presented in the previous 
paragraph. Instead of rotating the stimulus, we change its contrast as follows. Let us note Iq the stimulus 
peaked at and I n / 2 the one peaked at tt/2, the thalamic input to the hypercolumn of orientation takes 
the form 

/(i) = (l-V(t))io+-0(*)4/2 

where ip(t) is the function shown in figure 10. We check numerically, using the dynamics given by 
equations (7), that the hypercolumn stays in the cortical state TCq (see figure 11), despite the fact that 
the final stimulus corresponds to an orientation of 90 degrees. 

Discussion 

We now have a clear view of the functional impact of each parameter in the model. It turns out that 
the combination of mathematical and biological constrains basically fixes their values. We first notice 
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Figure 9. Plot of the position xq (shown in blue) of the peak of the stimulus / and the phase 
coordinate (p (shown in black) of the network, both as functions of time, on a logarithmic scale. Note 
that the stimulus drives the network into a state that is very close to its expected state when presented 
with a stimulus oriented at 90 degrees and that it stays there even after the input is switched back to a 
orientation stimulus. The parameters are the same as in figures 3 and 4. 
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Figure 10. Plot of the function ip allowing to vary the contrast of the thalamic input, see text. 
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Figure 11. The vertical direction represents the time, the horizontal one the orientation of the 
stimulus and of the network response. Left: Representation of the creation of the 90 degrees illusion. 
Middle: The stimulus starts as I to become a mixture of Jo and / 7r / 2 and finishes as a pure I^ii- Right: 
The response of the network shows that it stays in the cortical state TCq. 
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that the requirement for unimodal responses is a very strong constraint. Indeed, it implies that the first 
pitchfork bifurcation that occurs when A varies has to be the one corresponding to the first mode and 
this requires 

Ji > 0, J x > Ji Vi ^ 1 

Next, in order to actually see these tuning curves, the threshold 9 should not be too high: the condition 
shown in figure 2 is approximately J\ > 100 + 1 for £o — — 1. This in turn gives the range for the 
nonlinear gain A: it should be high enough in order for the model to produce tuning curves but smaller 
than the value for which the tuning curves saturate, hence do not vary with the input contrast anymore, 
in contradiction with the biological measurements. This means that A ~ -j-. The last relevant parameter 
is the width of the tuning curve. As shown in Supporting information, last paragraph, it can easily be 
estimated when 7T2 can be neglected which is often a reasonable assumption, see figure 6. It turns out 
to depend upon the ratio -?-. This closes the loop: we have set all the three parameters (A, J\, 9) and 

constrained the others (J, i >2). 

We now discuss the appearance of tuning curves (i.e. stationary solutions such that z\ ^ 0, \z\\ » 
\zi\, i > 1). What is the condition on the external input in order to produce such stationary solutions? 
Let us consider the case N = 1. We have seen (see figure 3) that if a tuning curve exists, then we have 
three stationary solutions and there is a value A c of the nonlinear gain A (roughly equal to 10 in figure 3) 
at which the two tuning curves disappear, we call it a turning point (it results from an opening of the 
pitchfork when e — 0). When varying e, we can look for the value of the nonlinear gain A (if there is 
one) at which a turning point occurs: it is an indication that tuning curves do exist for higher nonlinear 
gains. The TRILINOS package features the numerical continuation of the locus of the turning point and, 
starting with the turning point of 3, it can produce the locus of the turning points in the plane (e, A) as 
shown in figure 12. Above the blue curve, the stable response of the network is a tuning curve. Hence, if 

Locus of the turning points 
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Figure 12. Locus of the turning points as function of (e, A). Same parameters as in fig. 3. 



A > A c , even for no external input, e = 0, the network will be in a stable stationary state V*(x) which 
is a tuning curve whose tuning angle arg max x V* (x) is randomly selected. However, if A < A c until 
the contrast e has reached a certain value, the response will be untuned. It is more biologically relevant 
that the network operates in the regime A < A c , otherwise the neurons would have a high firing rate 
(around 60% of their maximum firing rate, see figure 3) even though no stimulus is present. Notice that 
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in figures 3, 4, we are in the case A > A c . Hence, for the parameters of figure 3, the working range of 
the nonlinear gain is A € [A m j n , A c ]. Notice that the condition on the threshold (see figure 2) gives the 
condition for A c to exist (in the limit e — > 0, the turning point converges to the pitchfork point). Also, 
elo lowers the threshold 8 and taking the differential of (13) w.r.t. 9, it can be shown that ^f e>0 < 
when eo < 0. Hence we expect that the local behavior around (0, A c ) in figure 12 to be quite general in 
the case £o < 0. This analysis hence provides a tighter constraint on the nonlinear gain A, it should be 
just below A c : A < A c . 

The fact that the cortical network shows two states corresponding to perpendicular orientations in 
response to a single stimulus can also be put in resonance with some published models of the cortical 
primary visual area (see [0] for a spatial network of Ring Models). Indeed, in this study of planforms 
in relation to visual hallucination, it may come as a surprise to the attentive reader that most of the 
planforms (in the cortical space) do not respect the good continuation principles of contours since adjacent 
hypercolumns show responses corresponding to orthogonal orientations. However once we agree that, for 
a hypercolumn, two orthogonal states are equivalent, this becomes perhaps less surprising. 

We relate our formalism to previous studies of recurrent models of orientation selectivity by first 
noting that the 90 degrees illusion was not reported in [ ] although they share the same assumptions as 
ours. 

In [15], the authors used a (voltage based) Ring Model in order to explain some of the features of the 
complicated spiking network of [2]. Although they used the non-saturating nonlinearity S(x) = max(a;, 0), 
they observed that narrowing the spatial extension of inhibition leads to mulimodal responses which they 
interpreted as neuronal illusions. This can be understood within our formalism: decreasing the spatial 
extent of inhibition introduces more Fourier terms (possibly with high values) in the connectivity J 
and can produce multimodal responses to a unimodal stimulus (see Finding the tuning curves, case 
N = 2). This type of nonlinearity (see also [ ]) cannot produce the 90 degrees illusion we have described, 
more generally, it is not possible to produce the 90 degrees illusions using non-saturating nonlinearities. 
Remember that the saturation arises when one takes into account the refractory period of neurons (see 
[16]), or more simply the fact that the firing rate of neurons is bounded. 

Under what conditions do the 90 degrees illusions survive in a network of Ring Models? Can we 
find similar illusions in more sophisticated networks and which experiments could confirm/invalidate 
our predictions? We just discussed the matter of a network of Ring Models with the study of [ ]. In 
[17], the authors used a generalization of the Ring Model with a very similar connectivity to explain 
the spontaneous activity observed in optical imaging recordings. They could not report the 90 degrees 
illusion because of their use of a non-saturating nonlinearity. If a detailed mathematical analysis of a 
modified version of their model to include a saturating rate function were conducted, it is very likely that 
the 90 degrees illusion would be predicted by the modified model. 

Finally, despite its ability to reproduce several experimental facts, the Ring Model lacks some anatom- 
ical data support because it does not use realistic cortical circuitry. Recently, M. Shelley et al. (see [18]) 
introduced a reduced system of a computationally intensive spiking neuron network model of a hyper- 
column with realistic cortical circuitry. Although they do not use a refractory period in their spiking 
model (hence, their reduced model has a nonsaturating nonlinearity), it could be interesting to look for 
neuronal illusions predicted by their model using the techniques developed in [13]. 

Conclusion 

We have pushed further the study, started in [13], of the mathematical properties of the Ring Model 
of orientation tuning and of some of their biological implications. This was achieved by taking into 
consideration the rich symmetries of the network. 
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For the first time to our knowledge in the field of neural networks, we have introduced the Orbit Space 
Reduction technique to deal with translation invariant connectivity kernels. This allowed us to find a 
suitable change of coordinates in order to remove the redundancy introduced by the symmetries. This is 
a generic technique that can be applied to many other problems in neuroscience. Using this reduction, 
we have shown that the exact shape of the connectivity function did not matter much as long as the first 
mode COS2 was the first to bifurcate. 

Our numerical continuation scheme has allowed us to discover another tuning curve encoded in the 
network that represents an orientation that is orthogonal to that of the LGN input. This neural illusion 
can be thought of as a ghost of the first pitchfork bifurcation that occurs when the sigmoid is centered 
(taking the value at the origin) and opens up when the bias on the sigmoid is removed. 

We have shown that it was possible to drive a hypercolumn to the illusory state by adding some 
dynamics to the LGN input: this gave rise to two dynamical illusions, one relying on rotating the 
stimulus, the other relying on changing its contrast. This is a strong prediction of the model that could 
possibly be tested experimentally even though this seems difficult given the fact that the Ring Model 
does not take into account the lateral spatial connectivity that is present in the visual cortex and allows 
different hypercolumns of orientation to interact with each other, but see the above discussion. 

It would be interesting to see if and how the illusions are modified when adding lateral spatial con- 
nections in a spatially organized network of Ring Models. 

Finally our approach leads to a near complete understanding of the role of each parameter in the Ring 
Model: the shape of the connectivity function through the weights Ji, the threshold 6, and the nonlinear 
gain A. 

Supporting information 

Numerical computation of the invariant functions 

Let us say a few words about the practical computation of the invariant functions Bo,a,b,c,d. As we 
are interested in the tuning curves, using the estimates of [13], we obtain that the L 2 -norm V' of 
the tuning curve is upperbounded by 14 for the connectivity function shown in figure 5. The relation 
II V II 2 = n ( v o + l^ 71 "! + if 71 ^) yields the estimation : 

,,, , ril 1-9-1 / (11 (91 Cauchy-Schwarz 3 

VI < N + y/JiQvPl + ki 2) |) + y/m(.\4 ] \ + \4 ] \) < -7= \\V f L < 24 
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for the same values of the parameters. The next step is to approximate the sigmoid S by a polynomial 
P on some interval [— a, a] where the value of a is chosen so that A V-H < a. As we need to observe 

L ' J II Moo — 

the first two pitchfork bifurcations, reached for the values Ai < A2 of A, see remark , we need at least 
A = A2 and, being a little bit conservative, we end up computing the solutions for A € [0,0.6]. This in 
turn requires a ~ 14. Note that the more accurate the approximation of S, the higher the degree of 
P with the consequence that some numerical instabilities may develop since this implies raising small 
numbers to high powers. 

P is then expressed in the basis of the Chebychev polynomials as P = J2i a iTi- The reason for this is 
that the Chebychev polynomials having rational coefficients, we can use, for example, the Groebner basis 
package of the symbolic computation package Maple to express the invariants Bq, B\, B^ as functions of 
(vq,tt). For example, we have: 



Maple 
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Bo = f f S[\v„ + \Vi(VJlz 1 e- 2 P t y + Vj;z 2 e- i P i y)]dy 

~ 2 

« f J P [\v + A3? (,/J[zie- 2piv + VJiZ2e- ipiy )] dy 

= f E 4 ", jf Ti [Xv + m (VT lZl e- 2 P*y + VT lZ2 e-^y)} dy 

~ 2 

fLiOifiiVo,*) 

"2" 

The computation of T^Uq, 7?) = / T, [Ai>o + A5R (-\/Ji-2 ; ie~ 2pii ' + \/~J~iZ2e~ 4piy )~\ dy is done automatically 

~~ 2" 

by the Groebner basis package but requires that the coefficients of the polynomial Ti be rational, not 
real. This justifies the Chebychev approximation. For a = 14 and an approximation error of 0.01 
(H5 — PHoo r 14 14 i < 0.01), it gives a polynomials P of degree 19. One important advantage of this 
method is that it does not require the vector (vq, 7?) to be on the Orbit Space to do the computations, i.e. 
we can compute B$, a, b 1 c, d even for values of (vo, 7?) that make no sense, e.g. ir± < 0, and then project 
the result on the Orbit Space. Note that the method is coherent since the results shown in figures 5 and 
G obtained by numerical continuation do satisfy A \\Vf\\ < 14, that is, are consistent with the numerical 

J J II II oo ' ' 

approximation. 

Numerical computation of the solutions of the nonlinear equations in the case 

N = 2 

In order to solve the nonlinear equations (14) for the tuning curves, we apply the strategy of [13]. The 
idea is to use a homotopy to solve the problem: we introduce a new parameter /i which translates 5*, 

i.e. Sfr) d = So + fi(-9 + 5(0)) where 9 is the threshold. Thus 5 (0 ) = 5 and S {1) = S - 9. This way 
we change the nonlinearity in (6) in order to find the TCs analytically (notice that this translation only 
affects the first equation of (6)). Indeed, when the nonlinearity is the centered sigmoid So we obtain the 
trivial solution V? — and we can also compute the values of the nonlinear gain A where the pitchfork 
bifurcations occur. We can then numerically continue this trivial solution with respect to the parameters 
(A, /x) to find the solutions of the equations with the "correct" nonlinearity, namely S — 9. We then 
simply take a slice of the output of the continuation program for fj, = 1 and obtain the dependency of the 
solutions w.r.t. the nonlinear gain A. This approach, though numerically intensive, is very convenient 
because it automatically gives the bifurcated branches. It also allows to compute some non-connected 
branches of solutions. This strategy relies on the library TRILINOS, see the acknowledgements below. 

Equivalent condition for the threshold 

Remember that our goal is to find a region in the plane (9, J\) where there exists a pair (A > 0, Vq) such 
that : 

, m .f 4 = eoS(Xv f )-9 

{Eh \l = A5'(A^)f 
We work out the case £o = 1, the other one being very similar. As 5' = 5(1 — 5), the second equation 
(E.2) becomes : 1 = ^5(1 - 5) ""'"fi^ 1 ) AA[/(_i _ [/) WIie re U d = v ! + 9. This quadratic equation 
in U has real solutions if and only if A J\ > 8, and they are given by 



-1±</1 



u ± '• x ' h 



s 
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We still have to verify that (E.l) is satisfied for at least one of these solutions. This yields an equivalent 
condition to (E) but with 3 variables instead of 4. 

For example, for U+, we obtain the equation in (A, Ji, 0) : U+ = —S(\U+ — X6). Using brute force 
computation for A £ [0, 30], 9 £ [0, 1], J\ £ [0, 10], we check when this is possible thereby obtaining the 
graphs shown in figure 2. 

The width of the tuning curves 

Lemma 0.1 If f(a,b) = — r -, then the width at half height of the tuning curve is equal to 



f(Xv£ - 8)A<fi{Ji) 



Proof. By definition, we look for the angle tp such that 



S 



A (v f Q +\Jit{j 1 -0\ = 25 A (v f + ^n{j 1 cos 2 (ip) - 0] . Setting a = X(v^ - 6) and b = 
\\/ir{ji, it follows that cos(2<p) = f(a, b). Hence the half-width is given by if = \ cos -1 f(a, b). D 
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